% This is the function that determines the optimal contract parameters for a country with commitment parameter mu_h, Total GDP of totalgdp1
% and an oilproduction of oilprod1.

%% Set the punishment parameter
function [result fval] = optparam(muh_input, nonoilgdp_input, oilprod_input, delta_input)


%% Determine the input parameters into the Weibull distribution, that generate a mean and variance as specified above

global mean sigma meanp beta muh mul pmuh pmul nonoilgdp oilprod delta xmax expon Vaut x0

mean  = 3.7063;  % Parameters of Lognormal distribution to fit oilprice to year 2000 - 2009.
sigma = 0.527;

logn  = @(z) z.*lognpdf(z,mean,sigma);       
meanp = quad(logn,0,10000);

display(['Mean Oil Price:   ', num2str(meanp)])

%% Parameters of the contract considered

nonoilgdp  = nonoilgdp_input;   % Empirically at the moment we are proxying for oilprod with the amount of government revenue 
                                % from hydrocarbon activities. Ideally, we would like to have a look at 
                                                
eta        = 0.5;               % Has to be between 0 and 1; Larger Values --> More risk aversion
expon      = 1 - eta;
oilprod    = oilprod_input;
muh        = muh_input;         % Cost of exproproiation, high state
mul        = 0.5*muh;           % Cost of exproproiation, low state
pmuh       = 0.90;              % Probability of cost of expropriation being high
pmul       = 1-pmuh;
delta      = delta_input;
xmax       = [350, 150, 200, 1]; % Initial Starting Values for Simultaneous equation solver
beta       = 0.9;   

%% Calculate V_aut - The autarky value

autarkyfunction = @(p) ((1/expon)*(nonoilgdp + (1-delta)*p*oilprod).^expon).*lognpdf(p,mean,sigma);
Vaut            = (1/(1-beta))*quad(autarkyfunction,0,10000);

display(['Vautarky:        ', num2str(Vaut)])
display(['------------------------------------------'])

%% Calculate Full-Insurance Utility

Util_Full_Insurance = (1 / (1-beta)) * (1 / expon) * (nonoilgdp_input + oilprod_input*meanp)^(expon);

%% Run fmincon over the parameters alpha and gamma

% Define Bounds
lb  =  [0        0];
ub  =  [meanp    1];
 
options1       = optimset('Display','on','TolFun',1e-6,'TolX',1e-5,'MaxFunEvals',10000);
[output fval]  = fminsearchbnd(@alphagamma, x0,lb, ub,options1);

alpha_f = output(1);
gamma_f = output(2);

% Now get resulting valus for p*, p**, Profit and Expropriation probability
SEValues = alphapstar1(alpha_f, gamma_f, meanp, beta, nonoilgdp, oilprod);

p1star_f    = SEValues(1);
p2star_f    = SEValues(2);
expprofitl  = @(z) ((1-gamma_f)*z -alpha_f + gamma_f*meanp).*lognpdf(z,mean,sigma);   % This is the expected profit to the company. Oil-price - Payments to country, which are given by: alpha + gamma(p - meanp)
ExpProfit_f = quad(expprofitl,0,p1star_f) + pmuh*quad(expprofitl,p1star_f,p2star_f); 
Prob_Exprob = pmuh * (1 - quad(@(z)lognpdf(z,mean,sigma),0,p2star_f)) + (1-pmuh) * (1 - quad(@(z)lognpdf(z,mean,sigma),0,p1star_f));

result = [alpha_f gamma_f p1star_f p2star_f ExpProfit_f Prob_Exprob]

display(['Alpha:          ',num2str(alpha_f,7)])
display(['Gamma:          ',num2str(gamma_f,7)])
display(['Prob. Exprop.:  ',num2str(Prob_Exprob,7)])
